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Abstract 

We have developed a general method for finding apparent horizons in 3D 
numerical relativity. Instead of solving for the partial differential equation 
describing the location of the apparent horizons, we expand the closed 2D 
surfaces in terms of symmetric trace-free tensors and solve for the expan- 
sion coefficients using a minimization procedure. Our method is applied to a 
number of different spacetimes, including numerically constructed spacetimes 
containing highly distorted axisymmetric black holes in spherical coordinates, 
and 3D rotating, and colliding black holes in Cartesian coordinates. 
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I. INTRODUCTION 



Black holes are among the most fascinating predictions in the theory of General Rela- 
tivity. The black holes most likely to be observed by future gravitational wave observatories 
(LIGO and VIRGO [|TJ) are those in highly dynamical evolutions, such as two colliding black 
holes. Moreover, events which are important for observations (i.e., events that occur more 
frequently and emit stronger radiation) are not expected to have a high degree of symmetry; 
for example, the inspiraling coalescence is a more probable scenario than the axisymmet- 
ric head-on collision of two black holes. The most powerful tool in studying such highly 
dynamical and intrinsically non-linear events is numerical treatment. 

The essential characteristics of a black hole are its horizons, in particular, the apparent 
horizon (AH) and the event horizon (EH). One needs to determine the location and the 
structure of the EH's in numerical studies to understand the properties of black holes, and 
indeed even to assert the existence of the holes. Algorithms for doing this have recently 
been developed |2]||. In contrast, the problem of determining the location of the AH in a 
general numerically constructed 3D spacetime has not yet been solved satisfactorily. The 
present paper represents a step in this direction. 

The apparent horizon is defined to be the outer-most marginally trapped surface [|J, a 
surface for which the divergence of the out-going null normal is zero (c.f., eq. (1) below). 
The surface is defined locally in time, in contrast to the EH, which can only be identified 
after the numerical simulation is complete. The AH, as a characteristic of black holes, can 
be used during the numerical construction of the spacetime. As discussed in a number 
of publications, apparent horizons are useful not only for studying the dynamics of black 
hole spacetimes ||, but also for use as an inner boundary in numerical evolutions of black 
holes ||-§|. The so-called apparent horizon boundary condition (AHBC) is currently being 
developed by many groups as a promising method for use in computing the long term 
evolution of 3D black hole systems. With AHBC one would like to be able to track the AH 
throughout the numerical evolution. For these reasons, it is important to develop efficient 
methods for locating apparent horizons in numerically constructed spacetimes. 

There are many well developed methods for determining the location of AH's in lower 
dimensional spacetimes [PHl3|l, e.g., in axisymmetry. The partial differential equation (PDE) 
defining a marginally trapped surface (eq. (1) below) reduces to an ordinary differential equa- 
tion (ODE) in the axisymmetric case, and the symmetry conditions also provide boundary 
conditions for starting the integration of the ODE. This simplifies tremendously the prob- 
lem, and enables the construction of efficient methods for finding the AH. However, as these 
methods rely strongly on the symmetry assumptions, they are not generalizable to 3D; going 
from 2D to 3D does not amount to simply adding one more spatial dimension. For the gen- 
eral 3D case, there is no symmetry and the AH surface to be determined is a closed surface 
(hence no boundary conditions for starting the integration) described by a non-linear ellip- 
tical PDE. At present there are no efficient algorithms for solving such a partial differential 
problem in general. 

We are aware of three independent efforts in determining the AH in the general 3D 
case. The first method is based on an expansion of the AH surface in terms of spherical 
harmonics, with the expansion coefficients determined by an integral equation. The equa- 
tion is then solved iteratively |L4l,|l5| . The second method attempts to solve directly the 
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elliptic PDE [|T^,|T^]. The third method [18| is based on expanding the closed surface in 
terms of orthogonal functions, in particular symmetric trace-free (STF) tensors, and using 
a minimization procedure to determine the expansion coefficients. Variations on this idea 
were explored by Brill and Lindquist jL9| and by Eppley |[20|| . The essential difference of 
this method with the first method is in the way the expansion coefficients are determined, 
and the use of STF tensor expansions for 3D Cartesian codes . The major advantage of this 
method is that the numerical solution of the minimization problem is much better under- 





we demonstrated that our method 



stood than the corresponding PDE problem. In Ref. 
can be efficiently applied to testbeds made up of analytically given data sets representing 
time-symmetric slices of spacetime. In Ref. the convergence of the symmetric trace-free 
tensor expansion for similar testbeds is studied in detail. In this paper, we follow up on our 
earlier work and give a more detailed discussion of this method. We also push the applica- 
tion of it in two directions: (1) Application of the method to arbitrary data sets which are 
not time-symmetric (non-zero extrinsic curvature), and (2) application of the method to 
data obtained in actual numerical evolutions of dynamical spacetimes. The spacetimes we 
studied include Schwarzschild and Kerr black holes, black hole plus Brill wave, and Misner 
two black hole spacetimes in full 3D. We demonstrate that our method is in principle appli- 
cable to general 3+1 spacetimes. Another method presently under development also uses a 
series expansion as in the present method, but it evolves the trial surface to the AH in an 
algorithm that combines elements of the AH finder of [14 and of a curvature flow method 
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In §H we discuss the formulation of our method, and the numerical algorithm in detail. 
Section [Til] gives the results of various testbed calculations, with §111 A| concentrating on 
initial data sets, and §111 B on spacetime evolutions. Section IV concludes with a discussion 



of where we stand in the construction of a general method for finding apparent horizons in 
3D numerical relativity. 



II. FORMULATION OF THE AH FINDER 
A. Basic Equations 

Defining to be the outward-pointing spacelike unit normal of a two-sphere S em- 
bedded in a constant time slice £ with unit normal n M , we can construct the outgoing null 
normal to any point on S as fc M = n M + s M . The surface S is called a marginally trapped sur- 
face (MTS) if the divergence of the outgoing null vectors vanishes V M & M = 0, or equivalently 



6 = DiS* + K ij s i s j - K = 0, (1) 

where G is the expansion of the outgoing rays, Di the covariant derivative with respect to 
the 3-metric 7^ , the extrinsic curvature of E, and K the trace of Kij. The AH is defined 
as the outer-most trapped surface. 

First suppose we are searching for the AH of a single black hole in spherical coordinates. 
As the AH is topologically a 2-sphere its position can be represented as 



F(M,r) = r 2 -/(M) = 0. (2) 
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The unit normal s l 



7 



^djF^dkFdeFy 1 ^, (3) 



can then be expressed in terms of the function f(8, <fi). Substituting (fj) and (|3|) into equation 
(H), one gets an elliptic equation for f(9, <p). Instead of solving this elliptic equation directly, 
we proceed by expanding / in terms of the usual spherical harmonics Y im : 

/(M) = E E F im Y^{9A). (4) 

1=0 m=-i 

Equation (|I]) then gives G in terms of the expansion coefficients F . The AH can be 
determined by finding the set of coefficients F £m which make = 0. 

For 3D codes in Cartesian coordinates, instead of the spherical harmonic expansion (f|), 
we choose to expand the trial surfaces in terms of symmetric trace-free (STF) tensors: 

F(x, y,z) = j2 (V - xlf - £ T Kl N Kl = 0, (5) 

i=l 1=0 

where x % are Cartesian coordinates and x the Cartesian points interior to the F = 
surface representing the horizon center. In expanding the function / in equation (|5|), we 
have adopted the notation from reference |f25f| . The STF tensors Tx t are coordinate in- 



dependent coefficients, the subscript Kt> is abbreviated notation for the vector product 
Nx e = nk 1 nk 2 ...rik n and the rii are unit directional vectors 

m = , =. (6) 

\>:] -..if 

To determine the set of which makes (5) the AH, we use a minimization procedure. 
In general, the AH surface is the outer-most surface represented by the set of and x 
which satisfy 

Y,W a Ql(^x )=0. (7) 

2 is semi-positive definite. W a is a positive definite weighting function which can be chosen 
to improve the accuracy of the minimization procedure depending on the construction of the 
trial surface. In this paper we do not investigate the use of this parameter and set W a = 1. 
With this we have converted the elliptic condition (1) to a standard minimization problem. 
There are efficient minimization algorithms to search the space of {Tk^Xq) coefficients for 
the surface closest to the apparent horizon among all test surfaces so parameterized. The 
strength of our method lies in that the minimization problem is much better understood 
than the numerical solution of the corresponding differential equation (1). 

An obvious potential difficulty of this method is that there is no guarantee that the 
summation converges, or converges rapidly. For black holes not in highly dynamical 
situations, we do expect the AH surface to be smooth and one needs only the first few 
lower rank tensors to find the surface accurately (as demonstrated in section [Hi]). This is 
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usually the case through most of the simulations, and our method will be more efficient at 
these times. However, we also expect that, for the cases we would like to simulate, there 
are often periods of time when the black hole (or holes) goes through highly dynamical 
evolutions, e.g., around the moment of the coalescence of two black holes. In such cases, one 
needs higher order STF expansions. Although STF tensor expressions are available in the 
literature, we find it convenient to include in our code a routine for the automatic generation 
of STF tensors to an arbitrarily high order. The procedure takes advantage of the fact that 
we can associate the with essentially x — Xq, y — yo or z — z in three dimensions. We can 
therefore construct all symmetric and independent permutations of 

N Kt =A(x~ x ) a (y - y f(z - z )\ (8) 

subject to the constraint a + (3 + 7 = £, where £ is the rank of the corresponding symmetric 
tensor, or equivalently the order of the multipole expansion, and A\ is the normalization 
factor making the RHS dimensionless. There are (£ + 1)(£ + 2)/2 such independent combina- 
tions. The combinations constructed in this fashion can be supplemented with the £(£ — l)/2 
independent conditions imposed on the symmetric permutations to make the rank £ tensor 
trace-free by contracting on any two indices. The function / can be expanded as 

f{x,y,z) = J2 E C k <{x-x o r k {y-y f«(z-z o y»>, (9) 

1=0 k 



where Ck are coordinate independent coefficients, and <> denotes STF combinations |25| . 
Partial derivatives diF and didjF needed to evaluate on the trial surfaces are then easily 
computed from (0) and ([5]). The simplicity of the form ([]) also allows one to easily construct 
the multipole expansion to take advantage of any symmetries present in the problem. For 
instance, our current implementation, by setting a flag, can invoke either the even or odd 
multipoles independently of the other, enforce axisymmetry, fix the surface centers x , or 
allow the most general parametric expansion. 



B. Numerical Algorithms 



The numerical problem is to find a set of parameters (Ck, x ) that minimizes the LHS 
of (0). Minimization techniques (such as conjugate gradient or quasi- Newton methods) that 
evaluate both the function and all its various partial derivatives are often preferred, as a 
means to increase the convergence rate, to those that do not require derivative information. 
However, because of the complexity of equation ([!]), in this first generation 3D AH finder 
we have chosen to implement a multi-dimensional method that does not require knowledge 
of derivatives, namely a direction set or Powell's method |[26|| . The method is based on 
successive line minimizations, whereby the function J2 ® 2 is successively minimized along 
different vector directions using the one-dimensional Brent's method with parabolic inter- 
polation. We find Powell's method to be generally robust, with a good convergence rate and 
computational speed, for surface functions parameterized by fewer than about fifteen or so 
parameters (see section |ITT|). Although solutions can still be found for parametrizations of 
higher order, the computational cost becomes overly excessive, especially when compared to 
the evolution cycle time. 
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A two-grid system is utilized to evaluate the expansion function O. The first grid (which 
we refer to as grid A) is the Cartesian-based computational grid on which the Einstein evolu- 
tion and constraint equations are solved, and the metric and extrinsic curvature components 
are defined. A second grid B is used to evaluate the surface function F(9, 0, r). This second 
grid need not be structured in the same way as the first. 

Assuming a guess position for the horizon center x l = (xq, yo, zo) on the structured 
Cartesian grid A, and some prescription for distributing points throughout this mesh to 
construct nodes for the second grid B, we can evaluate the expansion O on the nodes of 
grid B once f(6, 0) is determined. Here we choose to distribute the nodes uniformly over a 
sphere, centered on x l : the nodes are evenly spaced in both the polar and azimuthal angles 
to cover the full sphere < 9 < tt and < < 2n (or a single octant for the axi- and 
equatorial symmetric spacetimes). However, this procedure can easily be generalized to, 
for example, weight the node distribution according to the coordinate surface curvature, as 
might be desirable for highly distorted horizons. Along the radial direction, the nodes of grid 
B are placed uniformly with an inter-node spacing 5r typically equal to the cell resolution of 
grid A on which the Einstein equations are solved. We have also implemented a procedure to 
constrain the range of radii over which the solver searches for the AH, as might be desirable 
in cases where multiple trapped surfaces exist in the data set. In these cases, the radial grid 
spacing is set by 5r = {r max — r min )/N r , where N r is the number of radial nodes, and r min 
and r max are the lower and upper bounds of allowable radii. Representing the number of 
nodes on grid B as Ng x x N r , we set Ng = = 5 along the angular directions in a 
single octant, and N r = N along the radial direction, where iV is the minimum number of 
cells among the three orthogonal axes in grid A. 

Once the spherical grid B is constructed and centered on x l , the function F in equation 
(0) is evaluated on the nodes of grid B for a fixed set of coefficients (Ck,x l ) that the 
Powell routines compute. Along each radial line, a search is made for F = 0, by scanning 
from large to small radii, until the condition F(8,<f),r) x F(8,<j),r + 5r) < is met, and 
then linearly interpolating between adjacent neighbors to find the approximate Cartesian 
coordinates of the surface corresponding to F — 0. The extrinsic curvature and the metric 
functions and their derivatives are then evaluated at these positions by interpolating (either 
linearly or quadratically) from the computational grid A on which they are defined and 
used in equation ([[J) to evaluate 6 2 on the surface. The process of constructing a spherical 
grid B centered on x l , evaluating F on grid B, searching for the surface coordinates for 
which F = 0, interpolating the geometric data to the surface, and evaluating the expansion 
(H) on the surface is repeated through Powell's procedure until a minimum of @ 2 m the 
parameter space (Ck,x l ) is realized. 

To reduce the computational time spent in finding apparent horizons and to make the 
finder usable in real time evolutions of black hole spacetimes, it is important to implement 
the solver in a parallel fashion. Fortunately, the calculations performed at each node on 
the trial 2D surfaces are independent of the other node calculations. A natural parallel 
implementation, which we have adopted, thus distributes the different surface calculations 
to different computational processors, achieving a speedup of Ng x N^ compared to a purely 
scalar code. In addition, a large percentage of the cpu time is spent interpolating the field 
variables from grid A onto the 2D surface. To speedup this bottleneck procedure, we have 
written and implemented generalized algorithms designed to interpolate (both linearly and 
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quadratically) 3D data onto different nodes within a 2D surface in parallel. The computa- 
tionally intensive elements of the horizon finder have been optimized for both the Cray C90 
and the massively parallel Thinking Machines CM5 computers. 



III. CODE TESTS 

In order to test the basic solver described above, we have developed both 2D and 3D AH 
finders based on this minimization method. 

For the 2D finder, the AH determined can be directly compared to those obtained with 
the standard integration method. As testbeds, we used data obtained from a code developed 



by Bernstein et al. | 27fl . This code evolves a black hole distorted by an axisymmetric dis- 



tribution of gravitational waves (Brill waves) ||28|| . The black holes can be highly distorted 
by the incoming waves, leading to AH's with extremely prolate or oblate geometric shapes. 
In some cases the ratio of polar to equatorial circumference can exceed 10 2 . When these 
systems are evolved, the horizons undergo dynamic oscillations, eventually settling down to 
a Schwarzschild black hole at late times For such dynamic spacetimes, we compared 
the results obtained with our new AH finder algorithm with those from the AH finder con- 
structed using the standard ODE method ([T^,[J. For test surfaces with 16 coefficients, and 
using spherical harmonics YJ m as basis functions, we find that both methods produce the 
same results to within the given accuracy of the PDE solvers. 

As this paper is on a 3D implementation of these ideas for finding AH's, in the following 
we concentrate on results for the 3D case only. We have written a Fortran routine that 
implements the above method for a general 3D spacetime in Cartesian coordinates, and 
tested it on various spacetimes of interest. We discuss results from this code applied to 
various initial data sets containing one or more black holes (§ [111 A| ), and to evolutions of 



some of these black hole spacetimes carried out with our 3D "G" code ( §111 B|) . For all 
the 3D tests, we use the symmetric trace-free (STF) tensors as basis functions defined on a 
unit sphere. 



A. Finding Horizons in Initial Data Sets 

1. Schwarzschild Black Hole 

The simplest, most basic, test for any apparent horizon finder is the static Schwarzschild 
initial data. The 3-metric can be written in Cartesian coordinates as 

di 2 = (l + ^£) 4 (dx 2 + dy 2 + dz 2 ) , (10) 

where M is the mass of the black hole. The apparent horizon in this case is spherically 
symmetric and located at r = M/2. Although only the L = term should contribute to the 
surface that defines the apparent horizon, we tested the solver with a more general multipole 
expansion with the m/0 terms up to and including L = 6, a total of 28 coefficients. 

The computed surface is plotted in Fig. [1] for four separate cases with various black 
hole masses. The grid resolution used in each case was Ax = Ay = Az = 0.075M, using 
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25 3 cells. As expected, the surface is mostly defined by the i = contribution: The other 
higher order terms are small in comparison, roughly a factor of 10~ 8 smaller. We find typical 
errors in the horizon radius of order 0.04%, and the numerical surfaces are indistinguishable 
from the analytic solution in Fig. [3]. In each of the cases shown here, the finder converged 
to the correct surface in approximately 30 iterations. However, the number of iterations 
decreases significantly if fewer parameters are varied. For example, only 3 iterations are 
needed to converge if only the monopole term is varied, reducing the computational time 
by a factor of 60. On average, the cpu time scales approximately as N p , where N p is the 
number of parameters. Hence, the method becomes rather cumbersome for highly distorted 
horizons which can only be described with a high order multipole expansion. We discuss 
this important issue further in the following more elaborate tests. 



2. Misner Data 

The Misner initial data set represents two equal mass black holes initially at rest, and is 
defined by the 3-metric 

df = ^ [dx 2 + dy 2 + dz 2 ) , (11) 

where 

^ = 1 + E-i^(t L + — )> ( 12 ) 

^ smh(n/x) \ + r n r n J 

and 

± r n = \j x 2 + y 2 + [z ± coth(nyu)] 2 . (13) 

The parameter \i specifies the proper separation between the two holes and the total ADM 
mass of the spacetime. Apparent horizons in these data sets may consist of either a single 
surface surrounding both black holes if they are sufficiently close to one another (/i < 1.36) 
pOf , or two separate horizons located at the throats of the holes. 

In the cases where the two holes form a single encompassing horizon, the surface can 
be distorted significantly. To find the distorted horizons accurately, then, we need to keep 
higher order multipole terms, but a priori it is not clear how many terms will be required. 
In Fig. [2], we show the results of systematically increasing the number of axisymmetric 
expansion terms for the case \x = 1.2. In each calculation we use a 64 3 grid with Ax = 0.1 
and run the code on the 128 node partition of the CM5. We also show the result obtained 
with our 2D, axisymmetric code described in Ref. PI , which implements an independent 



ODE integration method ||. It is obvious that a high order expansion, up to L = 6, is 
required to accurately describe this surface which has a major to minor axis ratio of 1.5. We 
expect that even higher multipole expansions would be required for more distorted surfaces. 
In Table |, we show the number of iterations, cpu time and J2 © 2 as a measure of convergence 
for each of the expansion orders. We note that timings reported here and throughout this 
paper refer to the Thinking Machines CM5. 

For fi > 1.36, there are two separate and spherical trapped surfaces on the initial slice, 
centered off the origin at z — ±coth/i. In the cases we have tested (/i = 2.0 and 2.2), the 
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solver is able to locate both the center and radius of the offset horizons to an accuracy of 
better than 0.06% using general expansions to any order, L = to L = 6. We note that the 
iteration count (and hence cpu time) can increase by factors of between 2 to 10, depending 
on the total number of parameters, as compared to the cases in which the throat center 
coordinates are not allowed to vary. 



3. Black Hole plus Brill Wave 



The Misner initial data family just discussed provides examples of both single perturbed 
horizons and two separate but spherical surfaces with offset centers for testing. However, 
black hole horizons in highly dynamic spacetimes can be extremely distorted geometrically, 
and the horizon finder must be able to locate these as well. The black hole plus Brill wave 
initial data set is yet another solution that has been studied extensively in axisymmetry, 
and thus provides a useful testbed for highly distorted holes in three dimensions. This data 
describes the superposition of a black hole and a "doughnut" shaped Brill wave surrounding 
the hole. In spherical coordinates, the 3-metric takes the form 



2 • 2 

+ r sin 



(14) 



where q and ip are functions of r and 9 only. The function q is specified analytically as 
free data, and the Hamiltonian constraint is solved for the conformal factor ip. The initial 
extrinsic curvature vanishes due to time symmetry. 

This data set has been studied with 2D, axisymmetric codes [ 13 , 25 ] using a logarithmic 
radial coordinate 77 = ln(2r/mo), where mo is a scale parameter. In this coordinate system 
(77, 9, 0), the form of the q-i unction is written as 



Q 



a sin 



where we set n — 2, and 







+ exp 


fV~Vo\ 2 


exp 















(15) 



(16) 



We solve the Hamiltonian constraint for the conformal factor in our 2D code and then in- 
terpolate these solutions onto a 64 3 Cartesian grid with Ax = 0.15 to generate 3D data 
sets. The 3D horizon finder is tested against an independent solver developed for 2D calcu- 
lations 



m 



In Fig. H] we show the coordinate location of the apparent horizon for various parameters 
of the g-function. The sequence shown corresponds to different values of the Brill wave 
amplitude a, for fixed "shape" parameters having the values (r] ,a, n) = (0,1,2). Results 
from both the 2D and 3D calculations are shown. In the 3D case, we allowed searches up 
through L = 4, but for this case we include only the axisymmetric terms. In all cases shown, 
the finder was able to locate the correct surface to within a third of a grid zone in just 6 
iterations. 

Fig. ^ is a resolution convergence study of the case a = —0.75 in which the cell size is 
varied from Ax = 0.6, 0.3 and 0.15 with 16 3 , 32 3 and 64 3 grids respectively. The solver 
clearly converges to the correct surface quadratically with cell size. 
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For negative values of the Brill wave amplitude parameter a, the horizon is found off the 
throat, but for positive values below certain limits (depending on the shape parameters) the 
horizon is located on the throat at r = m /2. We are also able to locate the AH at the same 
level of accuracy for these cases. 

Although the surfaces shown in Fig. |3| appear to be almost spherical in coordinate space, 
the presence of strong gravitational waves can severely distort the horizons geometrically, 
much more so than the Misner data solutions described above. In Fig. ^ we show the 
(axisymmetric) geometric embedding of the two cases a = ±1, using the method described 
in Ref. ||. For all negative values of the amplitude parameter, the horizon is oblate. The 
a = +1 case is highly distorted geometrically into a prolate shape, with a ratio of polar to 
equatorial circumference of C p /C e = 4.28. 



4- Kerr Black Hole 

The calculations presented so far have tested the ability of the solver to find single or 
multiple horizons of spherical and highly distorted black holes in 3D Cartesian coordinates. 
However, the data in all these cases are time symmetric. The Kerr initial data set describing 
a rotating black hole has nontrivial extrinsic curvature, and thus provides another important 
testbed with a known analytic solution. The 3-metric for this spacetime in Boyer-Lindquist 
coordinates is given by 



2 (r 2 + a 2 Y - Aa 2 sin 2 9) sin 2 9 

dl 2 = P -dr 2 + p 2 d9 2 + ^ >- # 2 , (17) 

A p l 



where 



p 2 = r 2 + a 2 cos 2 9, (18) 

A = r 2 - 2Mr + a 2 . (19) 

In these coordinates, the non-vanishing components of the extrinsic curvature are 

K r( p = aM [2r 2 (r 2 + a 2 ) + p 2 (r 2 - a 2 )] sin 2 9/ (rp 4 ) , (20) 

K H = -2a 3 Mr cos 9 sin 3 9p~\ (21) 

where Kij = ip 2 Kij, and 

r = £±*ti***£l. (22) 

To construct this data in 3D, we first transform to an isotropic radial coordinate through 
the transformation 

/ M + a\ / M -a\ 

r = f 1 + — — 1 + — — , (23) 



2f J \ 2f 



as described in ||32|| . In this coordinate system, the coordinate singularities at r = M ± 



a/ M 2 — a 2 disappear. We then convert the metric and extrinsic curvature from isotropic 
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coordinates to Cartesian coordinates. The apparent horizon for the Kerr data is a coordinate 
sphere, located at 



VM 2 - a 2 

f = . (24) 

In Table [II] we show results from runs with various values of the rotation parameter 
a/M and expansion order L. These runs were all done with general expansions, without 
restricting to axisymmetry. For the tests performed here, Ax, Ay and Az were chosen to 
be approximately one tenth of the analytic AH radius, and the number of cells (32 3 ) was 
kept the same in all cases. Only the I = contributions are shown. As expected, the 
contributions from higher-^ terms are small, ranging from 10~ 9 for a/M = to 3 x 10~ 2 for 
a/M = 0.9. The reason that these terms, in the highly rotating not as small as 

in the Schwarzschild case is because the metric is not conformally flat, so the interpolations 
within the finder are not as accurate. Although we do not display embeddings of these 
horizons, we note that rotating black hole horizons can be extremely oblate, and for rotation 
parameters a/M > 0.866 a global embedding into Euclidean space is not possible. The 
horizon finder has no difficulty in locating these surfaces to within a grid zone, although we 
note, again, the high computational cost for large order multipole expansions as indicated 
by the a/M = 0.3 sequence of calculations. 



5. Transformed Schwarzschild Black Hole 

The above test cases, although treated in full 3D, have all been axisymmetric. In this 
section we test our finder on data which does not have any particular symmetry, and for 
which we can derive the correct solution in order to gauge the accuracy of the solver. To this 
end, we chose to find the apparent horizon on a Schwarzschild initial slice using a coordinate 
system in which the horizon surface does not appear to be axisymmetric. The coordinate 
transformation we use is 

/ l \-l/ 2 

f = rf (6, 0) = r (l + - sin 2 9 (cos 2 - sin 2 0) J (25) 

where r is the radial isotropic coordinate of (10). This / (6, 0) has surface modes of L = 2. 
The apparent horizon location is then defined by f = ^f{6,<p). 

The solver was allowed to search coefficients with terms up to L = 2 and L = 4. In all 
cases for the L = 2 tests, the finder successfully found the horizon to high accuracy. Each 
non-zero coefficient was accurate to better than a tenth of one percent, and the largest value 
found for a coefficient that was supposed to be zero was less than 2% of the smallest non- 
zero component. We note that better accuracy can be achieved if one uses full knowledge 
of the analytic form of the metric. When the finder was generalized to allow searches 
with coefficients up to L = 4, the horizon was again successfully found and to comparable 
accuracy. However with the large number of search parameters, some care was necessary in 
choosing the initial search direction of the Powell routine in order to avoid getting trapped 
at local minima in B 2 generated in the discretized representation of the spacetime. When 
the resolution of spacetime grid is increased, such local minima will diappear /change in 
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location, in sharp contrast to the actual G 2 = point (the actual location of the AH, which 
converges to a fixed location with increased resolution.) In Fig. || we show the apparent 
horizon found and the analytic solution in one quadrant of each of the coordinate planes. 
The numerical results shown in this figure were obtained with a run in which coefficients up 
to L = 4 were allowed to vary. 



B. Evolved 3D Data Sets 

The tests described thus far have shown that the horizon finder can locate horizons in 
a variety of distorted black hole spacetimes, but these have all been initial data sets with 
a somewhat restricted 3-metric and extrinsic curvature, and with a large part of the data 
prescribed analytically. In general 3D black hole evolutions, all metric and curvature com- 
ponents will be present, and the data will be contaminated with numerical inaccuracies 
generated during the course of evolution. The horizon finder must work under these condi- 
tions for it to be a useful tool in the numerical construction of spacetimes. In this section, 
we discuss results derived from the implementation of the horizon solver into our 3D "G" 



code 1 29] that solves the Einstein evolution equations in Cartesian coordinates. 



1. Schwarzschild 

The tests discussed here were carried out using multipole expansions up to and including 
the L = 4 terms. However, in order to save computational time, we restricted the search to 
axisymmetric surfaces. We have verified that a more general expansion does not change the 
results significantly. 

The results for a 3D Schwarzschild spacetime, evolved with both geodesic and maximal 
slicings with zero shift, are shown in Figs. [7| and [8] respectively. A 64 3 (130 3 ) grid with Ax = 
0.15 (0.2) was used for the geodesic (maximal) case with At = 0.25Ax. Only the i = 
contribution to the surface is plotted. The other parameters remain small (< 10~ 3 ) during 
the entire evolutions, as expected. In both cases, the surface locations are compared against 
the corresponding results from ID calculations. The two results agree to a small fraction 
of a grid zone throughout the evolution. The late time deviation in the maximal run is 
attributed to errors in the 3D spacetime evolution, which becomes inaccurate for t > 20M, 
although we note that the difference at the end of the run is still only about a grid zone. 

In both the geodesic and maximal cases, the AH finder accounted for a large portion of 
the total cpu time. For the geodesic case the code ran 160 timesteps, invoking the finder once 
every 0.3M in time (11 times total). In this case, the finder constituted approximately 90% 
of the cpu time. For the maximal slicing case the code was run for 1100 timesteps, calling 
the finder every 2M (14 times total). In this case, 40% of the cpu time was spent in finding 
the horizon. The horizon finder with a L = 4 axisymmetric multipole expansion is thus 
approximately 50-100 times slower than a single update cycle in the hyperbolic evolution 
(although we note that when the elliptic maximal slicing equation is solved, the relative 
performance in the solver improves significantly). A factor of roughly 10 can be gained by 
reducing the expansion order to L = 0. 
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2. Misner 2BH Collision 



A more difficult evolution scenario is to capture the horizon surfaces as two black holes 
collide and merge. Our "G" code is used to evolve the two black hole Misner data set, 



described in §111 A 2| , for a time sufficiently long that we can test the AH finder on this 
dynamic spacetime. 

To evolve the initial data, we use a zero shift vector and an algebraic lapse of the form 

a = a Q (1 + log 7) , (26) 

where 7 is the conformal 3-metric determinant and a a is the lapse on the initial time slice, 
which we take to be the Cadez lapse 



1 

a Q = - 



00 j 
^'^i ^ sinhn/x 



(27) 



The solution (|27j) solves the maximal equation in time symmetry with a = as a boundary 
condition on the throats. Algebraic lapses of the form ( |2"6| ) are singularity avoiding and 
produce evolutions similar to maximally sliced spacetimes [p9fl . 

The calculation is run on a 64 3 grid with Ax = 0.1 for a time of t = 10M using an 
expansion to order L = 4. As our interest is in testing the ability of the solver to locate 
the AH before and after the surface merger event, it suffices to evolve a data set with a 
low value for the Misner parameter, which is computationally less expensive in evolution. 
Here we show results for the fi — 1.5 case which has, as initial data, two coalesced black 
holes with a common event horizon, but two distinct trapped surfaces at the two throats. A 
common apparent horizon encircling both throats forms at time t ~ 1.6M. 

Figure || plots the surfaces found at each time the finder is called (t = 0, 2.5, 5, 7.5 and 
10M, where M is the ADM mass, and with the surfaces increasing in radius at the later 
times). At t = 0, the solver correctly finds the throat as the surface, centered off the origin. 
The finder can subsequently be prevented from locking onto the throat (which remains a 
trapped surface throughout the evolution) by restricting the range of radii over which the 
function F(9, (ft, r) is evaluated, and by resetting the center of the surface. Unfortunately, 
there is no analytic result to compare against the computed location of the AH at the 
later times. Instead, we overlay the 3D computed AH surfaces with the results from our 
2D axisymmetric code evolving the same initial data. The surfaces determined with these 
two different methods (3D minimization vs. ODE integration) applied to the two different 
constructions of the spacetime (3D Cartesian evolution vs. 2D axisymmetric evolution) 
basically agree with one another. Nevertheless, we note that in the figure the two sets of 
surfaces are not exactly the same in coordinate space, as they should not be, since we are 
using different lapse and shift functions in the two cases. However, the shift vector is of 
magnitude \(3\ ~ 10~ 3 in the region of the AH, so we do not expect this effect to exceed the 
grid spacing scale over the time interval of the calculation. Also, the algebraic lapse used in 
the 3D code is typically larger than the maximal slicing lapse in the 2D code by a fractional 
difference Aa ~ 0.02 near the AH. This should also have a small effect on the coordinate 
position of the AH except at late times. In fact, the maximum differences in the AH location 
between the 2D and 3D calculations is only slightly more than a grid cell. The differences 
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may also be attributed to inaccuracies in the evolution; as born out in computing the mass 
of the AH at different times. At t = 2.5M, after a common surface forms to surround both 
throats, we find Mah — 2.28, compared to the 2D result of 2.36. However, by t = 10M, the 
error in the mass of the computed surface increases to roughly 18% as the evolution becomes 
less accurate. These errors are attributable to numerical inaccuracies in the evolution and 
not the AH finder routine, in particular, they are unrelated to the multipole order used, 
since the surface becomes more spherical in time. 



3. Black Hole plus Brill Wave 

One final test of the solver in real time evolutions is a Schwarzschild black hole plus a 
Brill wave. The initial data is the same as that discussed in §111 A 3| with shape parameters 



(770, a, n) = (0, 1, 2), scale parameter m = 2 and amplitude a = —0.5. For this set of 
parameters, the ADM mass of the spacetime is M = 1.77m = 3.54. The Hamiltonian 
constraint is solved for the conformal factor in our 2D axisymmetric code, then interpolated 
onto a 66 3 Cartesian grid with Ax = 0.2 = 0.056M. The data is evolved with a timestep 
At = 0.25Ax = 0.014M, calling the horizon finder once every mo = 0.56M intervals of time. 
The shift vector is set to zero and the lapse function is computed from the maximal slicing 
condition. 



Figure [T0| shows the horizon shapes and locations in the 3D calculation at various times 
with intervals of 2mo = 1.15M, starting at t — 0. The multipole order used in this calculation 
is L = 4, with m = to reduce the computational time. The higher order expansion is 
needed here to describe the oblate shape of the horizons. For the most part, the solver 
required only about five iterations to converge when the previous solution is given as the 
initial guess to the finder. Also shown in Fig. |1(] are the corresponding surfaces found in our 
2D axisymmetric code. Once more we note that the slicing and shift conditions differ in the 
two cases, so we do not expect the surfaces to coincide precisely. However as the differences 
between the lapse functions and shift vectors are small (\(3\ ~ 10~ 3 and Aa ~ 0.07 in regions 
near the horizon except near the z— axis), the solutions do agree nicely. In and around the 
x — y plane, the solutions match to within half a grid cell. Greater differences, however, 
are found along the z-axis, where the two surfaces are displaced by a maximum of roughly 
two grid cells. This is attributed in part to a bigger Aa (~ 0.1 ) near the z-axis, which 
results from the imposed asymmetry in the lapse function due to the nearness of the outer 
boundaries, where we enforce the spherical Schwarzschild lapse as a boundary condition in 
the maximal equation used in the 2D simulation. 

Again, a more geometrical comparison or test of the solver is the mass of the surface 



found. The horizon mass is defined as Mah = y ^4ah/16vt, where Aah is the area of the 
surface. Figure [11] plots the AH mass as a function of time for the 2D and the L — 4, 3D 
evolutions. In both cases, the mass increases at first, as the gravitational waves fall into 
the black hole, reaching Mah ~ 0.997M at t ~ 6M. The masses in the two cases differ by 
only 0.1% at t ~ 6M. For comparison we also plot Mah for the surface found using a lower 
order L = 2 multipole expansion. At early times, when the horizon is most distorted, the 
L = 2 expansion is clearly not adequate to resolve the horizon shape, as evidenced by the 
AH mass which exceeds the ADM mass by about 1%. However, as the black hole settles 
down into a quasi-static state, the surface becomes more spherical and the L = 2 solution 
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approaches both the 2D and the L = 4, 3D results, differing from the 2D result by about 
0.35% at t ~ 6M. 



IV. CONCLUSIONS 

We have developed a promising general 3D method of finding the AH in a numerically 
constructed spacetime based on a minimization procedure. In this paper we have applied the 
method and demonstrated that it works for spacetime data which are not time symmetric, 
and data sets which are obtained in actual numerical evolutions. 

The major advantage of this method is that the minimization procedure is much better 
understood than the corresponding problem of solving the elliptic equation, and well-tested 
routines are available for solutions. The major drawback of the method is that, for AH 
surfaces which are not smooth, and which deviate significantly from sphericity in coordinate 
space, the higher dimensional minimization procedure can be computationally expensive. 
With our present implementation of the finder using Powell's routine, we are limited to 
searching for AH surfaces at every 50, or so, evolution cycles, instead of continually mon- 
itoring the AH throughout the numerical evolution, as would be our final goal. In future 
implementations of the method, we anticipate developing a more sophisticated minimization 
routine using derivative information to speed up the procedure, as well as other means to 
improve the robustness and efficiency of the code. Our code, which is optimized for both 
the C90 and massively parallel CM5 machines, together with documentation, will be made 
available on our servers |http:/ /jean-luc.ncsa.uiuc.edu| and |http:/ /wugrav. wustl.edu . 
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FIGURES 



FIG. 1. The coordinate location of the apparent horizon in the x-z plane for Schwarzschild 
initial data with various masses. The computed surfaces are indistinguishable from the analytic 
solutions, with deviations of order 0.04%. 

FIG. 2. The coordinate location of the apparent horizon for the fi = 1.2 Misner initial data. 
The solid line is the apparent horizon computed from our 2D code. The various broken lines are 
the surfaces obtained from the different multipole order expansions. The surface obtained from 
the L = 6 expansion is indistinguishable from the 2D result in this plot. We note that terms up to 
L = 6 must be included to find the AH accurately. 

FIG. 3. The coordinate location of the apparent horizon for 2D and 3D Brill wave plus black 
hole initial data with various values of the Brill wave amplitude. In all runs, the Brill wave is 
centered at rjo = and has a width of a = 1. Note that there is good agreement even though the 
initial data is not analytic and actually contains highly geometrically distorted black hole horizons. 

FIG. 4. A resolution study of the Brill wave plus black hole spacetime for the case (0,770,0") 
= (—0.75,0, 1). The surface converges quadratically with grid spacing to the correct location, as 
represented by an independent 2D calculation shown by the solid line. 

FIG. 5. The geometric embedding diagrams of the apparent horizons for the Brill wave plus 
black hole initial data with wave amplitudes a = ±1. The a = — 1 horizon is quite oblate, and 
deviates from sphericity even in its coordinate location. For the a = 1 case, the horizon is actually 
located on the throat, which is a coordinate sphere. However, the metric functions are highly 
nonspherical, leading to this very prolate geometry of the horizon surface. 

FIG. 6. Apparent horizon location in each of the coordinate planes for the transformed 
Schwarzschild initial data. The numerical data are represented by solid lines, and the analytic 
data by symbols. The surface is pure L = 2, although the finder was allowed to search through 
the L = 4 coefficients as well. The computed coefficients are found to roughly 0.1% accuracy. 

FIG. 7. The apparent horizon location for ID and 3D Schwarzschild evolutions using geodesic 
slicing. The ID data were obtained using 128 grid points and resolution Ar = 0.0375M, where 
M is the mass of the black hole. The 3D data were obtained using 64 3 grid points and resolution 
Ax = 0.075M. Only the t = contribution to the 3D apparent horizon is plotted but the other 
terms in the L = 4 series are negligible, as discussed in the text. 

FIG. 8. The apparent horizon location for ID and 3D Schwarzschild evolutions using maximal 
slicing. The ID data were obtained using 130 grid points and resolution Ar = 0.1M, the 3D data 
using 130 3 grid zones and resolution Ax = 0.1M. We note that the agreement is within a grid 
zone even to the end of the calculation, when the 3D evolution becomes inaccurate. 
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FIG. 9. Apparent horizon locations in the 3D evolution of the jjl = 1.5 Misner two black hole 
data. The surfaces are plotted at times t = 0, 2.5, 5, 7.5 and 10M, where M is the ADM mass 
of the spacetime, and the surfaces increase in radius with time. The evolution is performed using 
a 64 3 grid with resolution Ax = 0.1 and a multipole expansion of order L = 4. Also shown are 
the corresponding horizons found in the 2D axisymmetric evolutions. The surfaces from the two 
calculations agree nicely, although we note that an exact correspondence is not expected due to 
different shift and lapse functions. 

FIG. 10. Coordinate location of the AH found in the 2D and 3D evolutions of the Brill 
wave plus black hole spacetime. The surfaces are shown starting at t = with time intervals of 
2mo = 1.12M, where M = 3.54 is the ADM mass of the spacetime. Although we do not expect 
identical results due to different kinematic conditions in the two cases, the surfaces differ at most 
by little more than two grid cells. The evolution is performed on a 66 3 grid with Ax = 0.2. 

FIG. 11. Comparison of the apparent horizon masses computed from the 2D and 3D evolutions 
of the Brill wave plus black hole spacetime. The AH mass increases initially as the Brill wave falls 
into the black hole. By t ~ 6M, the mass approaches 0.997M and the 2D and L = 4, 3D 
results differ by just 0.1% at this time. We also show the corresponding masses computed from 
the surfaces found with an L = 2 multipole expansion. The low order expansion is clearly not 
adequate in resolving the surface at early times when the horizon is most distorted. 
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TABLES 



L iterations cpu time @ 2 



" " 2 ~ 1.2 min ' 8.3 x 10 2 

2 5 2.2 min 8.0 x 10~ 3 

4 9 4.6 min 7.9 x 10 -4 

6 16 12 min 1.0 x 10~ 4 



TABLE I. The effect of the (axisymmetric) multipole order L on finding the correct horizon 
surface is demonstrated here for the two black hole Misner data with \i = 1.2. The number of 
iterations and cpu time required by the solver is tabulated along with the expansion summed over 
all points on the surface. The timings were performed on the 128-node partition of the CM5, using 
a 3D grid of size 64 3 and a 5 x 5 mesh to cover the 2D horizon surface in a single quadrant. 



a/M 


L 


Ax 


r a 


r n 


Ar/Ax 


iterations 


cpu 





4 


0.1 


1.0 


1.0002 


0.003 


30 


2.4 hr 


0.3 


4 


0.1 


0.954 


1.0095 


0.556 


37 


1.0 hr 


0.3 


2 


0.1 


0.954 


0.9562 


0.022 


9 


6.7 min 


0.3 





0.1 


0.954 


0.9542 


0.002 


2 


0.5 min 


0.9 


4 


0.039 


0.436 


0.4748 


0.997 


27 


0.5 hr 



TABLE II. Performance measures of the horizon solver applied to the Kerr black hole data 
with various spin parameters a/M. Here L is the maximum multipole order, Ax is the 3D grid 
spacing, r a is the analytical position of the horizon, r n is the horizon location found by the solver, 
and Ar/Ax is the difference between the analytic and numerical locations normalized to the grid 
spacing. We also show the number of iterations and cpu time required by the solver running on a 
32-node partition of the CM5, using a grid with 32 3 cells and a 5 x 5 mesh for the surface. 



19 



Schwarzschild Initial Data 

I 1 1 1 I 1 1 1 I 



\ 



N ~ 



/ / 



\ 



\ 






X 



M=1 

M = 2 

M = 3 

M = 4 

_J i 

2 



Distorted Black Hole Convergence Study 




Schwarzschild Evolution, Geodesic Slicing 
2 o I i i i i i i i i i | i i i i i i i i i | i i i i 



.5 



1 .0 - 



0.5 

' D 
3D 



0.0 



12 3 

t/M ADM 



Schwarzschild Evolution, Maximal Slicing 



3.0 



2.5 



2.0 



1 .5 



1 .0 



0.5 



I I I I I I 


I 


I I 


I i i i i I 


i 


i i 


- 










- 
- 


— 










- 
- 




s ' 








- 


/. ■ ' 

- 

~ / ■'' 
■ — 

/ /' 






















/ '' 












- // 

- // 

I I I I I I 


I 






1 D 
3D - 

i i i i 


i i i 



5 10 15 20 25 

t/M ADM 



Misner Evolution, Data Evary 2.5M 




Distorted Black Hole Evolution 




